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Evidence for non-Fickian diffusion of a passive scalar is presented using direct simulations of 
homogeneous isotropic turbulence. The results compare favorably with an explicitly time-dependent 
closure model based on the tau approximation. In the numerical experiments three different cases 
are considered: (i) zero mean concentration with finite initial concentration flux, (ii) an initial top 
hat profile for the concentration, and (iii) an imposed background concentration gradient. All cases 
agree in the resulting relaxation time in the tau approximation relating the triple correlation to the 
concentration flux. The first order smoothing approximation is shown to be inapplicable. 
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I. INTRODUCTION 

In a turbulent flow the transport of a passive scalar 
is an important problem in atmospheric research, astro- 
physics, and combustion Passive scalar transport is 
also an important benchmark for more complicated tur- 
bulent transport processes such as turbulent magnetic 
diffusion and the alpha-effect in dynamo theory |3|, |J] , or 
turbulent viscosity and its nondiffusivc counterparts such 
as the AKA-effect H, and the Lambda effect 00. 

Modeling turbulent transport in terms of turbulent dif- 
fusion is known to have major deficiencies. For example 
turbulent transport is known to be anomalous, i.e. the 
width a of a localized patch of passive scalar concentra- 
tion may expand in time like tr 2 ~ t 13 , where (3 — 1 cor- 
responds to ordinary (Brownian) diffusion, (3 > 1 is su- 
perdiffusion, and (3 < 1 is subdiffusion |ij . Thermal con- 
vection, for example, has superdiffusive properties |10| . 
Turbulent transport is also known to have nonlocal and 
nondiffusive properties. One of the outcomes of this real- 
ization was the development of the transilient matrix ap- 
proach [Til Il2| which captures nonlocal transport prop- 
erties, although only in a diagnostic fashion In or- 
der to describe nonlocal aspects in a prognostic fashion, 
higher order spatial derivatives of the turbulent fluxes 
need to be included. These are best incorporated in terms 
of an integral kernel 01 ■ 

In the present work, however, instead of invoking 
higher order spatial derivatives, we follow the recent 
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proposal of Blackman and Field [Til UlSj to include an 
additional second order time derivative instead. This 
turns the diffusion equation into a damped wave equa- 
tion. Blackman and Field derived this equation from 
turbulent mean field theory by retaining triple correla- 
tions in the transport equation for the mean flux of a 
passive scalar. They assumed an isotropic turbulent flow 
and use a closure which relates triple correlations to dou- 
ble correlations [Trl IrH ITsl IT^ | . This approach is in some 
ways more elegant than the classical first order smoothing 
approximation 0, 0, |2(j , which breaks down because it 
assumes that the triple correlations are simply negligible. 
This approach also incorporates the momentum equation 
and, in magnetohydrodynamics, it therefore allows a nat- 
ural derivation of the feedback term of the alpha effect 
in magnetohydrodynamics 0, . 

Adding an extra time derivative in the equation for 
the turbulent transport of a passive scalar does certainly 
solve another long standing problem. Solutions to the 
diffusion equation are known to violate causality, be- 
cause the diffusion equation is elliptic and the propa- 
gation speed of a signal is infinite |2l|. This problem 
was originally discussed in the context of general relativ- 
ity 123, and more recently in the context of black hole 
accretion [2"3, 13 • The extra time derivative affects the 
modeling of turbulent transport most strongly at early 
times, just after having injected the passive scalar. This 
additional time derivative term tends to make the turbu- 
lent transport more ballistic at early times (correspond- 
ing to [3 ~ 2) . This property is well known in the context 
of standard Brownian motion. 

Non-Fickian diffusion has previously also been dis- 
cussed in various engineering applications, for example 
in diffusion problems in composite media I2M l26j| and 
in neutron transport problems in reactors |27j . which 
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are best modeled using non-Fickian diffusion. Here, 
a non-Fickian diffusion equation for particle transport 
arises by taking moments of the one dimensional Kramers 
equation, and approximating the second moment by the 
Maxwellian value 26, 28(. In these applications, however, 
turbulence is not considered. One exception is the recent 
work of Ghosal and Keller |2j| who derived a non-Fickian 
diffusion equation with the extra time derivative by going 
to the next higher order in an expansion of the under- 
lying integral equation. Comparing with data on smoke 
plumes in the atmosphere and on heat flow in a wind 
tunnel they find improved agreement with non-Fickian 
diffusion at small distances from the source. 

Given that the diffusion equation is now turned into a 
damped wave equation, one wonders whether oscillatory 
behavior is possible. Blackman and Field 0] find that 
oscillatory behavior is indeed present for long enough 
damping times but disappears for short damping times. 
For diffusion of a mean passive scalar, they argue that 
the oscillatory behavior is likely unphysical, and they 
use this to constrain their damping time to be of or- 
der of the eddy turnover time. However, the different 
numerical experiments presented below suggest that the 
damping time is about three times longer than the eddy 
turnover time. Furthermore, the simulations give direct 
evidence for mildly oscillatory behavior in a certain pa- 
rameter regime. 

The objective of the present paper is twofold. First we 
need to find out whether the existence of the proposed 
additional time derivative can actually be confirmed us- 
ing turbulence simulations. If so, we need to find out the 
magnitude of this extra term. Second, we need to study 
the range of modifications expected from this new term. 
In order to do this we consider numerical simulations of 
weakly compressible turbulence including the transport 
of a passive scalar. 

We begin by discussing the formalism that leads to 
the emergence of the additional time derivative in mean 
field theory. We then discuss the type of simulations car- 
ried out and present three numerical experiments that 
quantify the relative importance of the additional time 
derivative and that support the tau approximation for- 
malism. 



II. FIRST ORDER SMOOTHING VERSUS TAU 
APPROXIMATION 

A classic application of passive scalar transport is the 
diffusion of smoke in a turbulent atmosphere. If the 
smoke is injected in one point it will diffuse radially out- 
ward, so the mean concentration is expected to be a func- 
tion of radius r and time t. In that case it makes sense 
to consider averages over spherical shells, i.e. 

C(r,t) = — / C(r,8,<j),t)sm9d6d<j), (1) 
47r Jo Jo 



where C is the concentration per unit volume. Another 
application is the passive scalar diffusion between two 
parts of a slab that are initially separated by a mem- 
brane. In that case the mean concentration varies along 
the direction of the slab, say z, and then it makes sense 
to define horizontal averages, i.e. 

C(z,t) = — L- f * [ * C(x,y,z,t)dxdy. (2) 

J^x^y Jo JO 

This is also the type of average that is best suited for 
studies in cartesian geometry considered here. 

For clarity of the presentation here we ignore micro- 
scopic diffusion, in which case C satisfies the simple con- 
servation equation, 

W = -V.(IK7), (3) 

where U is the fluid velocity. The effects of finite micro- 
scopic diffusion will be discussed in the appendix. We 
now split U and C into mean and fluctuating parts, i.e. 

U = U + u, C = C + c, (4) 

and average Eq. ©, so we have 

flff 

— = -V-{UC + uc). (5) 

The challenge is now to find an expression for the concen- 
tration flux, uc = J- in terms of the mean concentration, 
C. The standard approach is to express the departure of 
the concentration from its average, c = C — C , in terms 
of its past evolution, i.e. 

c(x,t) = [ c{x,t')dt', (6) 
Jo 

where the dot denotes time differentiation and 

c= C-t= -V • (Uc + uC + uc-uc) (7) 

is the evolution equation for the passive scalar fluctuation 
obtained by subtracting Eq. (JSJ) from J21 • In the first or- 
der smoothing approximation or, which is the same, the 
quasilinear or second order correlation approximation Q , 
one ignores the terms that are nonlinear in the fluctua- 
tions, i.e. the terms uc — uc in Eq. J7J are simply omit- 
ted This is only justified if microscopic diffusion is 
large (but we have already assumed it to be negligible) 
or if the velocity is delta-correlated in time (which is also 
unrealistic). 

The terms that are nonlinear in the fluctuations would 
lead to triple correlations of the form mujdjC. Various 
authors have proposed to approximate triple correlations 
by quadratic correlations [l4llT(llT7lllMll3l which, in the 
present case, would be Uic/r; see Ref. [l5|. This is rem- 
iniscent of the Eddy-Damped Quasi-Normal Markovian 
approximation |30l l3l| , where fourth order correlations 
are approximated by third order correlations. This is 
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normally referred to as the tau approximation. In order 
to distinguish the two approaches, Blackman and Field 
[HI call the approach used in Refs [3 E IH E d the 
"minimal tau approximation" . In these approaches one 
calculates not J 7 , but instead its time derivative. In that 
case the time integration in Eq. Q disappears and one 
has 

d~T 



—— = u(x,t)c(x,t) + u(x,t)c(x,t). (8) 
dt 



This leads to the final equation 

dTi — Ti 

dt t 



(9) 



where r is some relaxation time and incompressibility 
has been assumed, i.e. djUj = 0. We shall now also 
assume isotropy, 

rms) where u rms is the rms 
velocity with u 2 ms = u 2 . The validity of Eq. © is clearly 
something that ought to be checked numerically using 
turbulence simulations. This is the main objective of the 
present paper. 

The other aspect is that the time derivative may not be 
ignorable in the final set of evolution equations. Thus, in 
contrast to ordinary Fickian diffusion, where the passive 
scalar flux T is assumed to be proportional to the mean 
negative concentration gradient (Fick's law), i.e. 



T = -K t VC (Fickian diffusion), 



(10) 



1 T U 2 
1 cor ""ri 



where Kt — 3 
fusivity and r cor 

in 



is the turbulent passive scalar dif- 
is some correlation time, one now has 



F=-K t VC-T— (non-Fickian), (11) 



where Kt = ^ru 2 ms . Equation l|lUfl can be reconciled 
only when time variations of the concentration flux have 
become small and if the correlation time r cor is identified 
with the damping time r. 

Applying dt + t _1 on both sides of 10, ignoring for 



simplicity a mean flow (U 
a damped wave equation, 

d 2 C 1 dC 



0), and inserting 1)11(1 yields 



dt 2 



dt 



(12) 



We note in passing that the extra term is in some ways 
analogous to the displacement current in the Maxwell 
equations. This is why this equation is also known in the 
literature as the Cattaneo-Maxwell equation [3^. The 
maximum signal speed is limited by u rms / V3- Assessing 
the importance of the extra time derivative is another 
objective of the present paper. 

The only ill-known free parameter in this theory is t, 
whose value is conveniently expressed in terms of the 
Strouhal number H, 



where kf is the forcing wavenumber or, more generally, 
the wavenumber of the scale of the energy carrying ed- 
dies. Here and elsewhere we consider it rms as a constant 
(independent of z and t). 

Some preliminary estimate of St can be made by con- 
sidering the late time behavior where Fickian diffusion 
holds. From Eq. (|10J) we expect that the decay rate of a 
large scale pattern with wavenumber k\ is 



A c = n t kl, 



(14) 



where Kt 



±tu 2 
3 v " 



is the turbulent diffusion coefficient. 



From forced turbulence simulations with initial mean 
flow or mean magnetic field patterns [33] , the decay rates 
of these patterns are well described by a turbulent kine- 
matic viscosity, vt, and a turbulent magnetic diffusion 
coefficient, 7y t , where both coefficients are about equally 
large with 

!/ t «?7 t «(0.8...0.9)xu rms /fe. (15) 

Applying the same value also to n t we obtain 

St ps (0.8. ..0.9) x 3 = 2.4.. .2.7. (16) 

This result is remarkable in view of the fact that in 
the classic first order smoothing approach to turbulent 
transport coefficients one has to assume St < 1; see 
Refs II i3. 

III. COMPARISON WITH SIMULATIONS 

In order to test the viability of the non-Fickian diffu- 
sion approach and to determine the value of St we have 
designed three different types of turbulence simulations. 
We first consider the problem of a finite initial flux, T , 
but with zero mean concentration, C = [l5j . Next we 
consider the evolution of an initial top hat profile for C 
and finally we investigate the case of an imposed uniform 
gradient of C which leads to the most direct determina- 
tion of t as a function of Reynolds number and forcing 
wavenumber. We begin with a brief description of the 
simulations carried out. 

A. Summary of the type of simulations 

We consider subsonic turbulence in an isothermal gas 
with constant sound speed c s in a periodic box of size 
2ir x 2ir x 2ir. The Navier-Stokes equation for the velocity 
U is written in the form 



DC/ 2 „, 

— = -c s 2 Vlnp + F v 

where p is the density, D/Dt 
advective derivative, 



/, (17) 
d/dt + U ■ V is the 



St — TM rms fcf, 



(13) 



v N 2 U + \ VV ■ U + 2S ■ Vlnp) 



(18) 
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is the viscous force where v = const is the kinematic 
viscosity, Sy = \(Uij + Uj,i) — \5ijUk,k is the traceless 
rate of strain tensor, and / is a random forcing function 
(see below). The continuity equation is 



dp 
di 



-V • (Up), 



(19) 



and the equation for the passive scalar concentration per 
unit volume, C, is 



dC_ 

dt 



-V • 



UC - pn c V [ - 
P 



(20) 



where kc — const is the diffusion coefficient for the pas- 
sive scalar concentration, which is related to v by the 
Schmidt number, 



Sc = v I kc- 



(21) 



Throughout this work we take Sc — 1. A nondimen- 
sional measure of v and hence kq is the Reynolds num- 
ber, which is here defined with respect to the inverse 
forcing wavenumber, 



Re = u rms / (vkf ) 



(22) 



The maximum possible value of Re depends on the reso- 
lution and the value of fcf. For fcf = 1.5 the typical value 
is approximately equal to the number of meshpoints in 
one direction. 

We adopt a forcing function / of the form 

f(x, t) = Re{Nf m exp[ik(t) ■ x + icf>(t)]}, (23) 

where x — (x, y, z) is the position vector, and — ix < 
4>(t) < 7r is a ((5-corrclatcd) random phase. The normal- 
ization factor is N — fo^kCg/St) 1 ^ 2 , with /o a nondi- 
mensional forcing amplitude, k = \k\, and St the length 
of the time step; we chose fo = 0.05 so that the maximum 
Mach number stays below about 0.5 (the rms Mach num- 
ber is close to 0.2 in all runs [I(|). The vector amplitude 
fk describes nonhelical transversal waves with \fk\ 2 = 1 
and 



/ fc = (fcxe)/Vfc 2 -(fe-e) 



(24) 



where e is an arbitrary unit vector. At each time step 
we select randomly one of many possible wave vectors 
in a finite range around the forcing wavenumber fcf (see 
below) . 

The equations are solved using the same method as in 
Ref. [37) . but here we employ a new cache and memory 
efficient code using MPI (Message Passing Interface) 
library calls for communication between processors. This 
allows us to run at a resolutions up to 1024 3 meshpoints 

MM. 



B. Finite initial flux experiment 

We consider first the example discussed by Blackman 
and Field In Fickian diffusion, if C = 0, there should 
be no flux, i.e. T — 0. Although this should in general be 
correct, one can imagine contrived situations where this 
is not the case, so it is an ideal problem to test whether 
the inclusion of the extra time derivative of the flux is 
at all correct and meaningful. Without this extra time 
derivative C would always stay zero. 

To explain in simple terms what happens, consider a 
situation where we have initially uniformly mixed white 
and black balls (so C = 0), but for some reason the balls 
are given an initial push such that the white balls move 
to the right part of the domain and all the black balls 
move to the left part of the domain. Then, after a short 
time, there should be a systematic segregation of white 
and black balls, in spite of continuous random forcing. 
Of course, this segregation survives only for a dynamical 
time, after which ordinary diffusion will begin to mix 
white and black balls. 

In order to set up such a situation in a turbulence 
simulation we assume that at t = the turbulence has 
already fully developed and then we initialize the passive 
scalar distribution according to 

C{x, y, z, 0) = Co smfciz. (25) 

Wrms 

Since u z = 0, and since the Reynolds rules Q are obeyed 
by our horizontal averages, we have C(z,0) = 0, but 
because v? z ^ 0, we have T z — v^c ^ 0. 

Numerically, we monitor the evolution of (C ) 1 / 2 , 
where angular brackets denote an average over z. This is 
to be compared with the analytic solution of the model 
equation (|12|l . Assuming that C(z,t) is proportional to 
exp(ifciz + Xt), the two eigenvalues are 



A±(fc 1 ) = -A ±AA(fci), 



(26) 



where 



U rms fcf . . 

Ao = ^sT' AA(fcl) 



A 



-U 2 

3 ""m 



M- (27) 



The solution that satisfies C(z, 0) = is 



(C 2 ) 1/2 = ^exp(-A t)sinh(AAt), 



(28) 



where A is an amplitude factor. Oscillatory solutions are 
possible (AA imaginary) either when St is large enough 
or, since St cannot be manipulated in a simulation, when 
fcf is small enough. According to Eq. I|lti|) we can estimate 

fcf/fci < 2 St/V3 « 3 (oscillatory behavior). (29) 

In the oscillatory case, AA is imaginary and so (C 2 ) 1 / 2 
is proportional to e~ A °*| sinwi|, where uj = ImAA. 

Note that the solution depends only on the combina- 
tion St/fcf, where fcf should be a known input parameter 



A 




Ur ms kft=G.6 Ur ms kft=l.8 ii rms fcff=3.3 
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FIG. 1: Passive scalar amplitude, (C } 1 ^ 2 /Co, versus time 
(normalized by u lnls kf) for two different values of fcf/fei. The 
simulations have 256 3 meshpoints. The results are compared 
with solutions to the non-Fickian diffusion model. 



TABLE I: Summary of fit parameters for the finite initial flux 
experiment. In all cases, the measured value of u Iuls — 0.23 
is used. Note that fcj flt ' is an independent fit parameter used 
instead of k{ to model the solution for a given value of k{. 
The range of wavenumbers used in the forcing function is also 
given. 





(range) 


fc f (flt) /fci 


st (fit) 


^(fit) 


1.5 


(1-2) 


1.0 


1.8 


0.21 


2.2 


(2...3) 


1.6 


1.8 


0.38 


5.1 


(4.5.. .5.5) 


3.8 


2.4 


0.18 



for a given simulation. However, in order to be able to fit 
the model to the simulation we have considered St and 
kf as independent fit parameters and refer then to the 
quantity k\ ^ . The results of our fits of the simulations 
to the models are shown in Fig.^ The corresponding fit 
parameters are listed in Table [I] 

We see that in all cases the Strouhal number does in- 
deed exceed unity. The resulting value is close to the 
value based on our simple estimate in Eq. 1)16(1. Second, 
oscillatory behavior of the solution is not only mathemat- 
ically possible for small values of fcf, see Eq. 1)26(1. but it 
is even physically realized in the solution for fef/fci = 1.5. 



C. Initial top hat function 

Next we consider the problem of an initial step func- 
tion. The advantage of such a profile as initial condition 
is that a broad spectrum of wavenumbers is excited. In 
order to avoid sharp jumps in the initial condition we 
choose a smoothed top hat function using the initial pro- 
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FIG. 2: C(x, 0, z) at three different times after reinitializing 
C according to Eq. I0DJ- fcf/fei = 1.5, Re LS = 400. 



file 



C(x, y,z,0) = h + h tanned 2 _ z 2 )], (30) 



where k z = 2 and d = 1 throughout this work. 

It is important to start the experiment at a time when 
the turbulence is fully developed. A visualization of C 
at three different times after reinitializing C is shown in 
Fig.H 

For Fickian diffusion the initial top hat function will 
broaden and develop eventually into a gaussian. As 
usual, for large enough values of the Strouhal number, 
wave-like behavior is possible and this would correspond 
to the initial bump splitting up into two bumps travel- 
ing in opposite directions. We have not been able to see 
this in our simulations so far. We have therefore decided 
to introduce as a quantitative measure of the departure 
from a gaussian profile the kurtosis, 



1 JCz 4 dz 
K ~V± JCdz ' 

where a quantifies the width of the profile with 



a 2 = 



JCz 2 dz 
JCdz ' 



(31) 



(32) 



For a gaussian profile we have k = 3, so we always plot 
k — 3. 

At early times, a 2 increases quadratically with t, but 
it soon approaches the linear regime, a 2 ~ t, until a 
saturates at a value comparable to the scale of the box; 
see Eq. <P|) . 

In Fig. El we compare the simulation results for a 2 and 
k — 3 with those obtained from the model ((12(1 using 
the same boundary conditions (periodic in z) and for the 
same values of u rms . For simplicity we solve Eq. 1(12(1 
numerically. However, similarly to the cases considered 
in § 1111 Bl we are unable to obtain good fits if we choose 
exactly the same values of k{ as in the simulation. There- 
fore, like in § 1111 Bl we treat kf as a fit parameter denoted 
by fc f (fit) ; see Table El 

There are characteristic departures in the values of a 2 
and k — 3 for the model compared with the simulations. 
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D. Imposed mean concentration gradient 




10 20 30 40 50 



model 
simulation 




10 20 30 40 50 



FIG. 3: Comparison of the evolution of a 2 and the kurtosis 
k — 3 for the non-Fickian diffusion model and the simulation. 
Note the good agreement at early and late times, but there 
are departures at intermediate times. The simulations have 
256 3 meshpoints. 



TABLE II: Summary of fit parameters for the initial top 
hat function experiment. In all cases, the measured value of 
Urms = 0.23 is used. Note that the values of St' flt ' are the 
same as those used in SII 



and the values of k* are now 



slightly closer to fcf than before. 



fcf/fcl 




st ( flt ) 


1.5 


1.3 


1.8 


2.2 


2.0 


1.8 


5.1 


4.6 


2.4 



Finally, we consider the case of a uniform gradient in 
the mean concentration. It is advantageous to split C 
into two contributions, 



C(x, y, z, t) = p(x, y, z, t)Gz + c(x, y, z, t), 



(33) 



where G — const is the imposed mean gradient of the 
concentration per unit mass (not unit volume) . Although 
C is now no longer periodic, this choice still preserves 
periodic boundary conditions for the departure c from the 
background profile pGz. Inserting Eq. I|33|) into Eq. I|20|) 
we have 



dc 



Uc - pn c V ( - ) - pn c Gz 



-pU z G, (34) 



where z is the unit vector in the z direction. The main 
advantage of this setup is the fact that we can now de- 
fine mean fields by averaging over the entire volume. We 
denote such averages by angular brackets. Note that 
(U) = 0, so U = u. The mean passive scalar flux is 
then (uc) and the triple correlation arising from (u z c) is 



Ti 



»V-(«c)). 



(35) 



Furthermore, there are triple correlation terms arising 
from the (ii z c) term via the momentum equation. The 
u ■ Vit term yields the triple correlation 

T 2 = ((uc) ■ Vu„), (36) 

and the pressure gradient term, Vp = Vlnp, yields 

T 3 = (cV zP ), (37) 

where p — In p can be regarded as a 'reduced' pressure 
and is related to the enthalpy. There is no correlation 
arising from the forcing term, because the forcing is delta- 
correlated in time. Furthermore, the contributions from 
the viscous and diffusive terms are small. Because of 
periodic boundary conditions, T\ + T 2 = 0, so the only 
contribution surviving from the sum of all three terms is 
T3. Thus, the final expression for r is 



(u z c) /{cV z p) 



(38) 



This could perhaps be explained by the fact that, espe- 
cially when k{ /k± is of order unity, the horizontal averages 
C obtained from the simulations are strongly 'contami- 
nated' by a small number of large eddies. Nevertheless, 
both at early and at late times the agreement between 
model and simulation is excellent. 

The results in § 1111 CI confirm our finding of § 1111 Bl that 
St is around 2 (or even larger). Again, this is large enough 
for oscillatory effects to appear when fcf jk\ is small. 



We note however that, on the average, the two con- 
tributions from the momentum equations cancel, i.e. 
T2 + T3 = 0. Therefore it is also possible to calculate 
t from the contributions of the passive scalar equation 
alone, i.e. 

r = (u z c) /(u z V ■ {uc)) . (39) 

We have calculated a series of simulations for differ- 
ent values of the Reynolds number as a function of k{. 
However, for a fixed value of v, and since k{ changes, 
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FIG. 4: Strouhal number as a function of kf/k\ for different 
values of ReLS- The resolution varies between 64 3 meshpoints 
(Re LS = 100) and 512 3 meshpoints (Re L s = 1000). 



FIG. 5: Visualizations of C on the periphery of the simu- 
lation domain at a time when the simulation has reached a 
statistically steady state, fcf/fci = 5.1, ReLS = 400. 



the Reynolds number, as denned by Eq. 122(1 . is not con- 
stant. Therefore we label here the curves by the value of 
the large scale Reynolds number that we define as 

Re L s = u ims /(vkx) . (40) 

The result is shown in Fig. 0] 

The resulting value of St depends weakly on kf and 
increases slowly with increasing kf. This dependence is 
weaker for smaller values of kf. As the Reynolds number 
increases, however, the range where St is approximately 
constant seems to increase. It is therefore conceivable 
that St converges to a universal constant whose value is 
around 3. 

Comparing with the work of Kleeorin et al. 0, 0] 
one has to note that the r approximation was originally 
formulated in fc-space (see also the early work of Orszag 
[3fj|). In Eq. 10, on the other hand, the r approxima- 
tion is applied directly in real space which may be the 
reason for minor differences. Nevertheless, under the as- 
sumption of Kolmogorov turbulence for k > kf, and no 
turbulence for k < kf , one finds that the Strouhal number 
is unity. Given that there can be further discrepancies 
arising from differences in the definition of St, we con- 
clude that their result is in broad agreement with ours. 

Since the simulations presented here are weakly com- 
pressible, comparison with incompressible theory may 
not be quite proper. If the assumption of incompress- 
ibility is relaxed, i.e. V • m ^ 0, there is an extra term, 
—UidjUj C on the right hand side of Eq. ©. In Eq. fl"^ 
this leads to an extra advection term, r -1 V • (t/ c ffC) on 
the left hand side. Here, U c s — U — ritV • u is a new 
effective advection velocity; see Refs. I n the sim- 

ulations presented here, the term uV • u is largest when 
kf/kf is small, but even then it is at most a few percent 




FIG. 6: Same as Fig.EJ but for h/ki = 1.5. 



of u^ ms fcf. This justifies a posteriori the neglect of com- 
pressibility effects in the interpretation of the numerical 
results. 

Visualizations of C on the periphery of the simulation 
domain are shown in Figures[S]and[n]for kf = 5.1 and 1.5, 
respectively [41| . Note the combination of large patches 
(scale ~ 1/fcf) together with thin filamentary structures. 
This is particularly clear in the case with kf/k\ = 1.5. 
The kinetic energy spectrum is close to fc~ 5 / 3 , but the 
passive scalar spectrum is clearly shallower (perhaps like 
k 1A ; see Fig.[7J|. These spectra are, as usual, integrated 
over shells of constant k = |fe| and normalized such that 
J °° E K (k)dk = i(u 2 ) and jg 00 E c (k)dk = \{c 2 ). 
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FIG. 7: Compensated kinetic energy and passive scalar spec- 
tra for the run with kf/ki = 2.2, ReLS = 1000. 

IV. CONCLUSIONS 

Two important results have emerged from the present 
investigation. First, the Strouhal number is generally 
above unity and may have a universal value between 2 
and 3 for forced turbulence. This implies that the clas- 
sical first order smoothing approach in invalid. Second, 
the triple correlations that are normally neglected are of 
comparable magnitude to the second order corrections 
that correspond to the passive scalar flux. The minimal 
tau approximation in which the two are assumed to be 
proportional to each other is shown to be justified. 

As was shown recently by Blackman and Field in the 
context of magnetohvdrodvnami cs |14l and then in the 
context of passive scalar diffusion [T]| , this leads to an ad- 
ditional time derivative in the mean field equation which 
then takes the form of a damped wave equation. Our 
work has now shown that when the forcing occurs on 
large enough scale (kf ^ 2fci) there is evidence for mildly 
oscillatory behavior. 

Among the various methods for determining the 
Strouhal number in a turbulence simulation, the ap- 
proach of imposing a uniform gradient of the passive 
scalar concentration is the most direct one in that no 
fitting procedure is needed. Using this approach requires 
however the firm knowledge that the functional form of 
the mean field equation is correct. This underlines the 
importance of the first two approaches where we were 
able to compare the evolution of various statistical quan- 
tities with those obtained by solving the model equation. 
The only shortcoming here is that we had to find not 
only the value of the Strouhal number, but we also had 
to allow kf to deviate (slightly) from the actual value 



of kf. Although the difference between the two is perhaps 
not unreasonable, one would like to have some theoretical 
understanding of this discrepancy. 

It is remarkable that in all three experiments the value 
of the Strouhal number depends only weakly on kf. This 
suggests that the relaxation time r decreases with in- 
creasing values of kf; see Eq. dl3fl . We also emphasize 
that St is similar in all three experiments, even though 
the wavenumber corresponding to the variation of the 
mean concentration changed a significantly. This sug- 
gests that r does not depend on the scale of the concen- 
tration, even though such a dependence is in principle 
being allowed for [l8L Il9ll30l | . 

The method used in the present paper to determine the 
Strouhal number from simulations can straightforwardly 
be applied to magnetohydrodynamics. In that case the 
magnetic field plays the role of the passive scalar gra- 
dient. Both satisfy very similar equations and in both 
cases a mean field can easily be applied while still retain- 
ing fully periodic boundary conditions. In both cases the 
closure approach of Blackman and Field predicts non- 
Fickian turbulent diffusion and hence the occurrence of 
an extra time derivative [T3 . [T5| . Their analytic approach 
and closure agrees reasonably well with our simulations. 
Another application would be to determine the role of an 
extra time derivative in connection with turbulent vis- 
cosity. In that case a mean gradient could be imposed 
using the shearing box approximation E3 • The first 
two methods described in the present paper should also 
still be applicable in that case. An obvious question that 
arises in this connection is whether non-Fickian diffu- 
sive properties could also play a role in attempts to find 
useful subgrid scale models for Large Eddy Simulations. 
The difficulty here is that one has to deal with averages 
that do not satisfy the Reynolds rules. Apart from this 
difficulty there should be no reason why an extra time 
derivative should not also be incorporated in such simu- 
lations. 
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